---
title: '[RC17]Cluster Annotation for Oligodendroglia'
author: "Nina-Lydia Kazakou"
date: "14/03/2022"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

# Set-up

### Load libraries

```{r message=FALSE, warning=FALSE}
library(SingleCellExperiment)
library(Seurat)
library(tidyverse)
library(Matrix)
library(scales)
library(cowplot)
library(RCurl)
library(scDblFinder)
library(Matrix)
library(ggplot2)
library(edgeR)
library(dplyr)
library(ggsci)  
library(here)
library(gtools)
```

### Colour Palette

```{r load_palette}
mypal <- pal_npg("nrc", alpha = 0.7)(10)
mypal2 <-pal_tron("legacy", alpha = 0.7)(7)
mypal3<-pal_lancet("lanonc", alpha = 0.7)(9)
mypal4<-pal_simpsons(palette = c("springfield"), alpha = 0.7)(16)
mypal5 <- pal_rickandmorty(palette = c("schwifty"), alpha = 0.7)(6)
mypal6 <- pal_futurama(palette = c("planetexpress"), alpha = 0.7)(5)
mypal7 <- pal_startrek(palette = c("uniform"), alpha = 0.7)(5)
mycoloursP<- c(mypal, mypal2, mypal3, mypal4)
```

### Load object

```{r}
oligos <- readRDS(here("Data", "Processed", "Oligodendroglia", "oligos_initial_srt.rds"))
```

## Plot marker genes for OPCs an oligodendrocytes
```{r}
DefaultAssay(oligos) <- "RNA"

Idents(oligos) <- "RNA_snn_res.0.05"
FeaturePlot(oligos, features = c("PDGFRA", "PCDH15", "MAG", "MBP"), label = TRUE)
Idents(oligos) <- "RNA_snn_res.0.3"
FeaturePlot(oligos, features = c("PDGFRA", "PCDH15", "MAG", "MBP"), label = TRUE)
Idents(oligos) <- "RNA_snn_res.0.5"
FeaturePlot(oligos, features = c("PDGFRA", "PCDH15", "MAG", "MBP"), label = TRUE)
Idents(oligos) <- "RNA_snn_res.0.7"
FeaturePlot(oligos, features = c("PDGFRA", "PCDH15", "MAG", "MBP"), label = TRUE)
```

The approximation models ran on this dataset point to the best (cleanest) resolution being RNA_snn_res.0.5 for the subseted oligodendroglia. However, based on the markers above we can see that in this resolution, cluster 5 (the one at the top), could be separated further, as half of it is PDGFRa+ (OPC marker) and the other half is positive for MPB+ (oligodendrocyte marker). So I have decided to use **RNA_snn_res.0.7** for my analysis as it seems to have more biological sense. 

```{r}
Idents(oligos) <- "RNA_snn_res.0.7"
```


# Cluster Annotation
```{r Cycling_progenitors, fig.height=6, fig.width=6}
FeaturePlot(oligos, 
            reduction = "umap", 
            features = c("MKI67", "TOP2A", "TPX2", "UBE2C", "CDK1", "PCLAF"), ncol = 2, label = TRUE)
```

```{r Oligodendrocyte_Astrocyte_Progenitor_Cells}
FeaturePlot(oligos, 
            reduction = "umap", 
            features = c("HOPX", "SPARCL1", "AQP4", "GFAP"), ncol = 2, label = TRUE)
```
OAPCs: similarly described in https://www.sciencedirect.com/science/article/pii/S2211124721001029?via%3Dihub

```{r Primitive-OPCs, fig.height=6, fig.width=6}
FeaturePlot(oligos, 
            reduction = "umap", 
            features = c("EGFR", "DLL3", "DLL1", "PPP1R14B", "PDGFRA"), ncol = 2, label = TRUE)
```

pri-OPCs: similarly described in https://www.nature.com/articles/s41467-021-20892-3


```{r OPCs, fig.height=6, fig.width=6}
FeaturePlot(oligos, 
            reduction = "umap", 
            features = c("PDGFRA", "PCDH15", "PTPRZ1", "COL20A1", "SOX6", "BCAN"), ncol = 2, label = TRUE)
```

```{r COPs}
FeaturePlot(oligos, 
            reduction = "umap", 
            features = c("SOX10", "MYRF", "NKX2-2", "MYC"), ncol = 2, label = TRUE)
```

```{r Oligodendrocytes,fig.height=6, fig.width=6}
FeaturePlot(oligos, 
            reduction = "umap", 
            features = c("MBP", "ZNF488", "PTGDS", "MAG", "MOG", "BCAS1"), ncol = 2, label = TRUE)
```


```{r message=FALSE, warning=FALSE}
oligos$ClusterID <- oligos$RNA_snn_res.0.7
current.ids <- c(0,1,2,3,4,5,6,7,8,9,10,11)
new.ids <- c("OPC_1", "pri-OPC", "OPC_2", "COP", "MAO", "OAPC", "OLIGO_1", "CyP_1", "Cyp_2", "OLIGO_2", "OPC_3", "SPARCL1+ OLIGO_3")
oligos@meta.data$ClusterID <- plyr::mapvalues(x = oligos@meta.data$ClusterID, from = current.ids, to = new.ids)
```

```{r}
oligos@meta.data$RNA_snn_res.0.01 <- NULL
oligos@meta.data$RNA_snn_res.0.05 <- NULL
oligos@meta.data$RNA_snn_res.0.1 <- NULL
oligos@meta.data$RNA_snn_res.0.2 <- NULL
oligos@meta.data$RNA_snn_res.0.3 <- NULL
oligos@meta.data$RNA_snn_res.0.4 <- NULL
oligos@meta.data$RNA_snn_res.0.5 <- NULL
oligos@meta.data$RNA_snn_res.0.6 <- NULL
oligos@meta.data$RNA_snn_res.0.8 <- NULL
oligos@meta.data$RNA_snn_res.0.9 <- NULL
oligos@meta.data$RNA_snn_res.1 <- NULL

oligos@meta.data$Dataset <- "hES-derived Monolayer Oligodendroglia"
oligos@meta.data$Author <- "Kazakou et al."
```


```{r}
head(oligos@meta.data, 2)
```

```{r}
DimPlot(oligos, group.by = "ClusterID", pt.size = 0.9, label = TRUE, cols = mycoloursP[6:40]) & NoLegend()
```

### Check how many cells from each timepoint make up each cluster for my chosen resolution.
```{r}
n_cells <- FetchData(oligos, 
                     vars = c("ClusterID", "orig.ident")) %>%
  dplyr::count(ClusterID, orig.ident) %>%
  tidyr::spread(ClusterID, n)
head(n_cells, 5)

write.csv(n_cells, here("outs", "Processed", "Oligodendroglia", "no.ofCellsperCluster.csv"))
```

### SOme more plots
```{r}
DimPlot(oligos, group.by = "Sample", pt.size = 0.9, label = TRUE, cols = mycoloursP[4:40]) & NoLegend()
```

```{r fig.height=5, fig.width=13}
DimPlot(oligos, group.by = "ClusterID", pt.size = 0.9, label = TRUE, cols = mycoloursP[4:40], split.by = "Sample") & NoLegend()
```

```{r}
DimPlot(oligos, group.by = "Phase", pt.size = 0.9, label = TRUE, cols = mycoloursP[15:40]) & NoLegend()
```

```{r fig.height=5, fig.width=13}
DimPlot(oligos, group.by = "ClusterID", pt.size = 0.9, label = TRUE, cols = mycoloursP[4:40], split.by = "Phase") & NoLegend()
```

```{r}
DimPlot(oligos, group.by = "RC17_AllCell_ClusterID", pt.size = 0.9, label = TRUE, cols = mycoloursP[1:40]) & NoLegend()
DimPlot(oligos, group.by = "RC17_AllCell_RNA_snn_res.0.8", pt.size = 0.9, label = TRUE, cols = mycoloursP[1:40]) & NoLegend()
```

```{r eval=FALSE}
saveRDS(oligos, here("Data", "Processed", "Oligodendroglia", "hES-derived_monolayerOL.rds"))
```

```{r}
sessionInfo()
```